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AN INVESTIGATION OF THE OPTIMUM DESIGN 


AND FLIGHT OF ROCKETS* 

By John D. Bird 
Langley Research Center 

SUMMARY 

An analysis using classical variational methods was made whereby significant 
physical characteristics of rockets were determined in an optimal sense. A two-point 
boundary -value problem was formulated which, when solved, provided the significant 
physical parameters and flight path for a two-stage rocket which would deliver a given 
payload to orbit with minimum total initial mass. Numerical solutions were carried out 
for nine cases by an iterative method. The results show how the stage masses and 
thrusts should be proportioned for various structural efficiencies in the fuel tanks and 
rocket motors. Analyses of the type carried out for this investigation should be useful 
in the sizing and, generally, in the design of rockets. 

INTRODUCTION 

Many applications of variational methods exist wherein the flight paths of rockets 
or aircraft are optimized. In general, a solution is found for the flight path that will 
enable the vehicle to deliver the maximum payload to orbit, arrive at the terminal condi- 
tions in minimum time, or accomplish the mission in some other way that is regarded as 
being most efficient. Another, and in some ways a more complete problem than is nor- 
mally treated by variational methods, involves solution for both the flight path and the 
vehicle configuration which will deliver a given payload to desired terminal conditions 
with the object of minimizing the initial mass of the vehicle or some other quantity of 
significance. 

The purpose of this paper is to present an investigation of this latter problem for 
a rocket application as completely as seems reasonable for a generalized analysis of 
relations between performance and structural design parameters and to illustrate the 
effects of substantial changes in structural efficiency as might be representative of 

*The information presented herein was included in a dissertation submitted in par- 
tial fulfillment of the requirements for the degree of Doctor of Science in Aerospace 
Engineering, University of Virginia, Charlottesville, Virginia, July 1967. 



various degrees of conservatism in design practice. Solutions were found for the most 
efficient two-stage rocket from the standpoint of minimum initial mass that would place 
a given payload in a circular orbit about the earth at an altitude of 370 400 meters 
(200 int. n. mi.) with a continuous burn of the motors. In these solutions the burn times 
and thrusts of the two stages were determined as well as the steering of the thrust vector. 

In order to formulate this problem, a brief analysis of engineering data on the con- 
struction of rockets was made to determine the significant parameters in defining the 
mass of a rocket. This analysis enabled expression of the mass of a two-stage rocket in 
terms of certain constants representative of good construction and of the burn times and 
the thrusts of the two stages. 

Several papers have been written that contribute to the solution of problems of this 
type. (See refs. 1 to 5.) Bliss, in reference 1, gives the basic multiplier rule of the cal- 
culus of variations. Denbow, in reference 2, derived necessary and sufficient conditions 
for a minimizing arc for problems having free corner points. His solution is restricted 
to the case where the state variables are continuous at the corners. Mason, Dickerson, 
and Smith (ref. 3) have extended the multiplier rule formulation of Denbow to the case 
where discontinuities at various points in the minimizing arc may be treated as functions 
of the state and independent variables at the ends and the points of discontinuity of the 
minimizing arc. This extension allows the optimization of the relative sizes of rocket 
stages. References 3, 4, and 5 give examples of this type of problem and its solution. 

Lockheed Missiles & Space Company has had good success with problems of this 
type by use of the method of steepest descents. Reference 6 gives examples of solutions 
for problems similar to those solved in this paper. With the exception of Lockheed's 
work, the thrust levels in the various stages and the structural elements between stages 
have not been included in analyses of this type. 

For this investigation the methods of the classical calculus of variations as given 
in reference 3 are used to derive the necessary conditions in the form of differential and 
algebraic equations that must be solved in order to obtain the desired result. The 
methods result in a two -point boundary -value problem wherein part of the conditions on 
the variables is known at each end of the problem. Solutions are presented for nine cases 
that are representative of a range of structural design practices from very heavy tanks, 
motors, and interstage structures to zero mass structural elements to show the influ- 
ence of design practice. A comparison is made with a nominal design for the purpose of 
illustrating the benefit of configuration optimization. This nominal design was chosen for 
the beginning of the iterative solution process. 
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SYMBOLS 

cross-sectional area of rocket, 37.12 meters^ 
drag coefficient, 0.5 
admissible solution 
force vector 

functions of given arguments 

gravitational acceleration at earth surface, 9.807 meters per second^ 

H-function of Bliss and Pontryagin 

altitude, r - r e , meters 

specific impulse of rocket fuel, 400 seconds 

unit vectors of polar coordinate system 

functional of given form 

rocket motor structural mass fraction 

rocket interstage structural mass fraction 

rocket fuel tank structural mass fraction 

constant multipliers 

mass, kilograms 

mass of payload, 4207 kilograms 

designates number of variables or functions in series, appendix B 
maximum dynamic pressure, newtons per meter^ 
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r e 

t 

*a 

*c 

*1 

*2 


*3 

*1 + 

V 

*2 + 


V 


u 


u 


a 


u 


X 1 


x 2 


x 3 

x 4 


radius vector 

magnitude of radius vector, | r |, meters 
radius of earth, 6 371 203 meters 
time, seconds 

independent variable, appendix B 

value of t at a corner, appendix B 

time at initiation of rocket flight 

time of first- stage burnout, staging time 

time of second-stage burnout at orbital injection 

refers to instant after t^ 

refers to instant prior to tg 

refers to instant after tg 

refers to instant prior to tg 

control force vector 

magnitude of control force vector, |u|, newtons 

unit control force vector 

magnitude of radius vector, r, meters 

time rate of change of magnitude of radius vector, r, meters per second 

angular displacement of radius vector, 9 , degrees 

time rate of change of angular displacement of radius vector, 6 
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x 5 


rocket mass, m 


Xg inclination of thrust vector to local horizontal, /3, degrees 

X 7 magnitude of thrust vector, u 

x 20’ x 21 ,x 22’ x 23 functions defined by equations (15) 
a angle of attack, degrees 

/3 inclination of thrust vector to local horizontal, xg, degrees 

y flight-path angle relative to local horizontal, degrees 

6 angular displacement of radius vector, xg, degrees 

Aj Lagrange multipliers, functions of time 

Aq constant multiplier 

v constant for atmospheric density exponential function, 1/6705 per meter 

p(h) mass density of atmosphere as a function of altitude 

p e mass density of atmosphere at earth surface, taken as 

1.752 kilograms per meter^ 

first-order differential function 

1 p algebraic function 

A dot over a symbol denotes a time derivative; the form F x . denotes a partial 
derivative; a repeated subscript indicates a summation over that index. 

The notation | j indicates the magnitude of a vector. 



ANALYSIS AND PROBLEM FORMULATION 


Assumptions 

In the solution of this problem a number of assumptions were made in an effort to 
reduce a rather complex problem to a form for which numerical solutions could be made 
with a minimum of programing difficulties. More complete problems are substantially 
more difficult in programing, checkout, and solution. The assumptions made for this 
investigation are summarized as follows for convenience. Some of these assumptions 
are discussed later. 

(1) The initial mass of the rocket is taken as the significant parameter for 
optimization. 

(2) Point mass equations of motion are used. 

(3) A planar trajectory is used. 

(4) The rocket is assumed to have a constant cross-sectional area. 

(5) The aerodynamic drag coefficient of the rocket is assumed to be independent of 
Mach number and angle of attack, and the lift-curve slope of the rocket is assumed to be 
zero. 

(6) The earth rotational velocity is assumed to be zero. 

(7) An exponential function of altitude is assumed to represent adequately the 
atmospheric density. 

(8) The effect of the atmospheric pressure on the rocket motor thrust is assumed 
to be zero. 

In this investigation the assumption is made that the initial mass of the rocket is a 
significant parameter for minimization, and the problem formulation and numerical solu- 
tions are carried out on this basis. To some extent this assumption is valid but actually 
the problem is more complex because of the fact that rocket motor development is a long 
and expensive process whereas fuel tank enlargement is generally a simpler process. 

As a result, when the payload capacity of a rocket must be increased the least expensive 
step is to increase the fuel capacity and move into a less efficient mode of operation inso- 
far as fuel consumption is concerned. Much additional thought needs to be given to the 
choice of an appropriate payoff quantity or quantity for optimization. A criterion based 
on cost would certainly seem to be most significant. 

Point mass equations of motion are employed in this analysis. The necessary atti- 
tude changes should be slow and of little consequence in the overall trajectory determina- 
tion. This assumption is common in trajectory analysis. Moreover for computational 
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simplicity only motion in the vertical plane is considered which reduces the motion to 
one in two degrees of freedom in a plane passing through the center of the earth. The 
equations of motion are derived in appendix A. The rocket is assumed to have a con- 
stant cross-sectional area throughout the numerical solution process. Physically, this 
assumption requires that the fuel tank size be changed in the solution for the optimum by 
varying the length of the tank rather than by changing its diameter. An exponential func- 
tional representation was employed for the atmospheric density. This representation 
depends on the altitude and a single density gradient constant and has been shown to be 
reasonably accurate. (See ref. 7.) 

Rocket Structural Mass 

An analysis of available literature on liquid fuel rockets was made in order to 
establish what is presently considered to be good rocket design practice and to formulate 
rocket mass as a function of the significant parameters of the problem. As a result of 
this analysis the following relationships were established: 

m rocket motor hardware = k M u 

m fuel tank = k T m fuel 

m inter stage structure = ^S m elements above J 

where kj^, krp, and kg are constants determined by design practice. The first rela- 
tionship indicates that the mass of all hardware directly associated with a liquid fuel 
rocket motor is directly proportional to the magnitude of the vacuum thrust for which the 
motor was designed. The second relationship is generally well known among rocket 
designers and indicates that the mass of all hardware directly associated with the con- 
tainment and transfer of liquid fuel or propellant is directly proportional to the mass of 
the fuel to be contained. This mass includes tanks, bulkheads, baffles, and major piping. 
The third relationship indicates that the mass of the interstage structure is directly pro- 
portional to the mass of all components mounted above it. This relationship indicates 
that the mass of a supporting structure is directly dependent on the mass to be supported 
as might well be assumed. 

An element of mass not associated with rocket motor thrust or fuel mass, but 
rather more mission dependent, is associated with each stage of a rocket. This element 
of mass is associated with mission control and telemetry, in addition to other factors. 

No effort was made to account for this mass in the analysis although it is significant. 


a) 
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These relationships are naturally only valid for rockets constructed and erected in 
the earth gravitational and atmospheric environment. The interstage and tank structure 
might be appreciably less if the rocket could be built and erected in a lunar environment. 
No dependence on rocket acceleration is shown in these relationships. This is undoubtedly 
associated with the rather substantial structural mass needed for satisfactory handling 
and erection of a rocket under the earth gravitational and atmospheric conditions and to 
the approximate nature of the analysis and the limited amount of data readily available 
for analysis. Exponential relationships seem to give a more exact representation of the 
structural weights, but the linear relationships used here were felt to be adequate. 



Figure 1.- Rocket components and masses at ignition time. 


By use of these three relationships 
(eqs. (1)) and the knowledge that the thrust of 
a rocket depends on the specific impulse of 
the fuel used and the time rate of fuel consump- 
tion according to 


u= -Ig e 


dm 

dt 


( 2 ) 


an expression can be derived for the mass of 
a two-stage liquid fuel rocket as a function of 
time and certain parameters such as the con- 
stant thrust of each stage, the burning time of 
each stage, the payload mass, and the mass 
fractions determined here. This relation 
between thrust and mass -flow rate does not 
include a correction for atmospheric pressure. 
The relationship between rocket mass and 
parameters is based on the two-stage rocket 
structure shown schematically in figure 1. In 
this figure the mass of all elements of the 
rocket is given at ignition time t^ in terms 
of the significant parameters of the problem. 
The general relationship for the mass of the 
rocket for any time during burn is 
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m(t) = (1 + k s) [o + k s) m P + t 1 + k T) Ugt ff ? (*3 - *2) + k M u g |- a g e 

+ (1 + k T) St ^ -(*2 " *1) + k M u stage lj “ S ^® 6 "(* " *1) (*1 

m (t) = (1 + k s )m p + (1 + k T ) ^i (t 3 - t 2 ) + k M u stage 2 - ^§pl(t - t 2 ) 




< t < U 


} ( 3 ) 


(t 2 < t < t 3 ) 


Thrust 


where m p is the payload mass, is the ignition time, t 2 is the staging time, and 
tg is the time of burnout. The thrust magnitudes u are labeled for the appropriate 
stage. This formulation for the mass of the rocket is used in the derivation of the neces- 
sary conditions for a minimizing arc. 

Formulation of the Two-Point Boundary -Value Problem 

The problem formulation of concern here involves the determination of the mini- 
mum mass, two-stage rocket capable of placing a specified payload in a circular orbit 
at an altitude of 370 400 meters. The essential 
parameters of the rocket such as thrusts and 
stage burn times are to be determined along 
with the trajectory in space. The rocket 
coordinate geometry is shown in figure 2. 

The necessary conditions that must be sat- 
isfied in order to assure an extremal solu- 
tion and, hopefully, a minimum solution to the 
variational problem at hand may be derived 
from the multiplier rule stated in appendix B. 

The differential constraints of the prob- 
lem may be put in first-order or state -vector 
form by appropriate substitutions. To 
accomplish this, the following terms are 
defined: 
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*1 = r(t) 
x 2 = r(t) 


x 3 = 0(t) 

X 4 = 0(t) 
X 5 = m(t) 


x 6 = j8(t) 


X7 = u 


where, also, Xg has been introduced as the direction of the thrust vector u as given 
in figure 2 and xrj has been introduced as the magnitude of the thrust. The magnitude 
of the thrust is fixed during any one stage; hence, its differential equation assumes a 
simple or null form, specifically X7 = 0. In line with the definitions. 


X 1 = x 2 


x 3 = x 4 


Utilizing these definitions together with the equations of motion derived in appen- 
dix A and the formulation developed for the rocket mass (eq. (3)) enables the development 
of the following differential constraint functions for use in the formulation of the neces- 
sary conditions. These functions correspond to cpj of appendix B. 


0j(t,X,x) = Xj - X£ = 0 

4> 2 (t,x,x) = x 2 - Xl x 4 2 - ^ sin x 6 + g e r e 2 ^-I-j + C D A | p(h) ^(x^ 2 + x 2 2 ^ = 0 
0 3 (t,x,x) = x 3 - x 4 = 0 

0 4 (t,x,x) = x 4 + - ^ cos x 6 + C D A ^( x 1 2x 4 2 + x 2 2 ) V2 = 0 


>(4) 


^ 5 (t,x,x) = X 5 + - — =0 

i§e 

= x^ = 0 
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Incterms of the state vector xj the rocket mass as derived previously becomes 

x ( t *) 

(1 + k g )m p + (1 + krp) 7 jg "^(*3 ' *2) + k M x 7 ( t 2 + ) 


x 5 = (i + k s ) 




t 1 + kr r)~ljr (*2 ■ 4 l) + k M X 7 ( t l + ) 


c 7( t l + ), t _ t x 

ig G P 'u 


(ti < t < t 2 ) 


> ( 5 ) 


x 5 


= (i + kg)m p + (1 + k T ) 7 jg^ ^(*3 - *2) + k M x 7 ( t 2 + ) ‘ 7 jg“'( t ' M 

(t 2 < t < tg) 

The function F defined in the multiplier rule of appendix B is given as 


F = Xfy 


= x l[*l -* 2 ] 


+ Ac 


x 2 - Xl x 4 2 - -1 sin X 6 + g e r e 2 -±^+0^1 p(h) ^(x^ 2 + x 2 2 


>i 


X 3 [^3 “ x 4] + X 4 |*4 + ” xpq cos x 6 + C D A £ p < h) £~( x l 2x 4 2 + X 2 2 ) J 

[V j^j +X 6[°] + X 7[ x 7] (6) 


+ X 


Then from Euler's equation (eq. (Bl)), 4r (^- ] = the following equations may be 

y ox^ J 


obtained for arcs of continuous x^: 


X 1 - -X 2 


H 2 + 2g e r e 2 ( ri) - C D A ^ ^( x 1 2x 4 2 + x 2 2 ) ^ x l x 4 

1/2 


2x„x 4 x 7 

n COS X R 

Y 2 y 2y ® 

X 1 X 1 x 5 


- c d a ?^ 4( Xi2x42 +X22 ) 

- C d A i p(h) + x 2 2 ) 1/2 x l x 4 2 - C D A I ^ sf ( x 1 2x 4 2 + x 2 2 ) 1/2 ] < 7a > 
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*2 - -XiOG - x 2 JjcnA | p(h) i-(x! 2 x 4 2 .f x 2 2 ) 1/2 - CjjA i p0>) ^(x^ 2 + x 2 2 ) 7 x 2 J 


^ 2 + x 2 V 1 / 2 x 

v 4 +x 2 J x 2 


X 3 = 0 


(7b) 


(7c) 


* 

X 4 = 


~ x 2 [^ x l x 4 “ C D A |pO») ^|( x 1 2x 4 2 + x 2 2 ) ^ x 1 2x 4 ] “ X 3 H 

C D A jP( h ) ^( x 1 2x 4 2 + x 2 2 ) V2 - C E ^( x 1 2x 4 2 + x 2 2 ) 1/2 x 1 2x 4 ] 


(7d) 


x 5 - -x 2 


_ _1 gin x 6 + CjjA | p(h) — |{x 1 2 x 4 2 + x 2 2 ) 
x 5 


+ CjjA |p(h) _|(x 1 2x 4 2 + x 2 2 ) 
x 5 


x< 

x 5 

1/2 




- X. 


x l x 5 


■ COS Xc 


sin x c 


X 7 " " X 2 X5 ' a 4 Xl x 5 + ig e 


cos Xg x 5 


(7e) 


(7f) 


The Euler equation associated with the component xg serves to define the thrust 
direction for an extremal solution in terms of two of the multipliers. This result is 


T x 7 1 T x 7 • 1 

0 = X« cos x R + sin x R 

2 L x 5 6 J 4 [?1 x 5 6 J 


which can be written 


X 2 

tanx 6 = -xj 


( 8 ) 


This result allows for two choices for xg. This ambiguity is resolved by the necessary 
condition on the Weierstrass E-function which is equivalent to requiring that the function 
H = XjXj be a maximum. 
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Also since H is explicitly independent of t and since from the multiplier rule 


dH _ 9H_ 
dt 8t 

H is a constant along the arcs of continuous Xj of the problem. 


( 9 ) 



For application of the end and corner conditions of the multiplier rule (eqs. (B2) 
and (B3)), the function G is written and contains the optimized quantity for this problem, 
the initial mass, and, in addition, the necessary end and corner constraints of the problem. 
The function G is given as 


G = X Q g + l ^ 

which, in terms of the problem, becomes 

G = X 0 x 5( t l + ) + 




( 10 ) 


where the quantity to be minimized g is the initial mass of the rocket and where the 
constraining relations are given as 

*1^1 = h | x s(V) - i 1 + k s) m P - k T ^"^3 - *2) - ^(V) 

HH = l 2 | X 5(V) - x 5 ( t 2 + ) - k S { 1 + k s) m P + { 1 + k T)-^-^( t 3 - *2) + k M x 7 ( t 2 + ) 


* 3^3 _ ^3 


*4^4 " ^4 


{*> ' °) 

( t l + ) " r initial^ 
= h ( x 2( t l + ) " initial] 
*6^6 = ^6 1^3 ( t l + ) " ^initial} 


(Hb) 

(He) 

(Hd) 

(He) 

(Ilf) 
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= Irj |^ X 4 (tj j - ^initialj 

(Hg) 

*8^8 = *8( x l(V) “^(V)) 

(llh) 

*9^9 = l 9 ( X 2(V) " ^(V)} 

(lli) 

*10^10 = *10 (* 2 ") ■ X 3( t 2 + )J 

(Hj) 

*11^11 = *11 (*2 ) ' x 4( t 2 + )j 

(Ilk) 

*12^12 = *12 ( x l(V)- r final] 

(HI) 

h3 ^13 = l 13 {^(V) ' r finalj 

(11m) 

*14^14 = l 14 f x 4(V) " ^finalj 

(Hn) 


The constraining relation ^ is an expression for the final mass of the rocket in terms 
of the defining parameters and is obtained from the expression for rocket mass as a 
function of time. The constraining relation is an expression for the mass disconti- 
nuity at the staging point t 2 and is also obtained from the expression for rocket mass. 
The constraining relations !//g, and if/rj express the known initial con- 
ditions of the problem. The relations and express the known final 

conditions of the problem. The relations i//g, \f/g, and i/'u express the fact 

that the components of the state vectors x^, Xg, Xg, and x^ are continuous at the 
staging point tg. 

The end and corner conditions given in equations (B2) and (B3) require that certain 
expressions be zero at the points tj + , tg", tg + , and tg", where the superscripts + 
and - indicate conditions immediately after and prior to the designated point in time. 
These expressions give relations between the variables of the problem at these points 
and enable the end conditions of the problem to be established. Figure 3 gives the 
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conditions associated with the present 
problem. The variables that are subject to 
solution in this problem are listed for the 
points tj+, t2~, t2 + , and tg". The com- 

ponent X7, the thrust magnitude, is free at 
all four of these points subject to the differ- 
ential constraint £7 = 0 because X7 is 
held constant on each continuous arc or roc- 
ket stage. The values of certain components 
are fixed at the launch and orbital condi- 
tions, of course. 

Substituting the values of F and G 
of this problem (eqs. (6) and (10)) into the 
end and corner relations (eqs. (B2) and 
(B3)) gives the following relations that must 
be satisfied. The multiplier Xq is taken 
to be 1 since the relations are homogeneous 
in the multipliers. For the point tj_, 


Orbit 



Figure 3.- Free variables at various points in the trajectory. 


*5 ‘ *2( l l + ) _ 0 

h - x s(*i + ) ' 0 

h • x 4( l i + ) = 0 

1 - x 5 (t,+) = 0 

h * H (‘i + ) - 0 


( 12 ) 


for the point t2, 


*8 + x l(*2~) = 0 

(13a) 

-*8 - x l(V) ’ 0 

(13b) 
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*9 + x 2 (*2 ) “ 0 
-l 9 - ^( t 2 + ) = 0 


Z 10 + X 3( t 2 ) = 0 
~*10 “ X 3 (* 2 + ) = 0 
l ll + X 4 ( t 2 _ ) = 0 
~ l \\ ‘ X 4^2 ) = 0 
l 2 + X 5( t 2 _ ) = 0 
“ X 5( t 2 + ) = 0 

Xrj(t 2 ) = 0 



(13c) 

(13d) 

(13e) 

(13f) 

(13g) 

(13h) 

(13i) 

(13j) 

(13k) 

( 131 ) 

(13m) 
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The relations involving Zg, Zg, Z^q, Z-q, and Z 2 and the values of Xj on each 
side of t 2 show that the multipliers Xj, X 2 , Xg, X 4 , and Xg are continuous at tg. 
The relations involving Z 4 , Z 5 , Z 0 , Z 7 , Z 3 , Z 12 , Z 13 , and Z 14 reflect the initial and 
final conditions on the state variables of the problem. 

Substituting to eliminate Z^ and Z 2 in the appropriate relations and defining 
four new quantities x 2 q, x 24 , x 22 , and x 23 give 



" x 7( t 2 + ) “ 0 


The end and corner conditions are summarized in brief form in table I, with the 
exception of the conditions on x 2 q, Xgp x 22 , and x 23 ; numerical values for the end 
conditions of the launch trajectory are given in table n. Table I shows the conditions 
that are fixed at given values and the conditions that are free to be determined by the 
problem solution. 

In integrating from tg“ to tj + it can be seen that x-^tg - ^, x 2 ^tg - ^, x 4 ^tg - ), 
Xg^tg X 7 ^t 2 and | are specified and that Xg^tg X^^tg )* 

X 4( t 3 _ )’ X 5( t 3")’ x 7( t i + )> x 7( t 2 + )» * 2 ~> ^ t 3~ are free to be chosen to satisfy 

conditions on x^*), x^t^), Xgftj+J, x 4 (tj + ), X 5 (t 1 ' t '), x 2Q , x 21 , x 22 , and 

x 2 g. The value of Xg is determined at all points in time by the specification of t 2 ~, 
tg _ , x 7 , and x 7 (t 2 + ). Note that the number of free conditions corresponds to the 
number of conditions to be satisfied. The quantities t 2 " and tg" occur in the equa- 
tions in the forms t 2 - tj and tg - t 2 which are equivalent in information content 
because tj^ is zero. 
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TABLE I.- END AND CORNER CONDITIONS ON TRAJECTORIES 


Variable 

Values of variable at specific times 

*i + 

t 2 _ 

t 2 + 

t 3 _ 

t 

0 

Free 

V 

Free 

x i 

r initial 

Free 

x i(V) 

r final 

x 2 

r initial 

Free 

x 2( t 2") 

^final 

x 3 

^initial 

Free 

x 3( t 2") 

Free 

x 4 

* x 5 

^initial 

Free 

rt- 

W 

1 

Rfinal 

x 7 

Free 

x 7( t l + ) 

Free 

x 7( t 2 + ) 

X 1 

Free 

Free 

X l(V) 

Free 

x 2 

Free 

Free 

^(*2 ) 

Free 

X 3 

Free 

Free 

x 3( t 2 _ ) 

0 

X 4 

Free 

Free 

X 4( t 2 _ ) 

Free 

X 5 

1.0 

Free 

X 5( t 2 _ ) 

Free 

x 7 

Free 

0 

Free 

0 


*xg(t) is determined explicitly by a choice of t 2 ~, 
V> x 7 (t 1 +), and x 7 (t 2 +). 

TABLE H.- NUMERICAL VALUES OF END CONDITIONS 


OF TRAJECTORIES 


x l(tl+) or r initial . . . 

6 371 203 m 

x 2 (ti + ) or r initial . . . 

0 

X 3( t l + ) or ^initial • • • 

0 

x 4( t l + ) or ^initial • • 

0 

x l( t 3 ) or r orbital * * 

6 741 712 m 

x 2 ( t 3 _ ) or ^orbital • * 

0 

X 3(V) or ^orbital • • 

Free 

MV) or Vbital • • 

0.00114023 rad/sec 


In the numerical solution of a problem, values of the quantities 

X 2(V)’ X 4(V)’ h{V)> *2 - h> *3 " * 2 ’ x 7( t l + )^ ^ x 7( t 2 + ) 

an effort to cause the quantities x^(tj + ), Xg^tj 4 '), x 4 ( t i + ) 


x 3( l 3")- x l( t 3 _ )* 
are chosen in 

’ X 5( t l + )» x 20’ 
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Xgj, — X 22 > and Xgg to achieve required values. These last quantities are obtained by 
integrating the differential equations of constraint and Euler equations in reverse and 
assembling the necessary quantities from the results. The quantity x i(t]+) must be 
equal to the radius of the launch point, X^(t^ + ) must be equal to 1.0, and the other 

and x 3 (tj + 

remainder of the solution is independent of these quantities. 

Method of Solution 

A differential correction procedure of simple form, the Newton-Raphson method, 
was used to obtain numerical solutions to the two-point boundary -value problem. In this 
procedure the differential equations that define the evolution of the state and multiplier 
vectors in time are integrated along a nominal path starting with some known and some 
assumed conditions at one end of the problem, and then along adjacent paths corresponding 
to perturbed values of each of the assumed parameters. The incremental changes in the 
boundary and other required conditions are obtained, and a matrix of partial derivatives 
is constructed from which a system of linear equations is obtained that relate changes in 
the boundary conditions to changes in the parameters free for variation. These equations 
are solved by matrix inversion, and corrections are made to the parameters free for 
variation so that a better match of required conditions results. Because of the nonlin- 
earity of the problem, a diagonal weighting matrix was employed to limit the degree to 
which the various boundary conditions were corrected in one iteration. The necessity 
for employing this weighting matrix to avoid instability of the process of iteration to the 
optimum solution was recognized by Robert E. Smith, of the Langley Research Center, 
who carried out the programing for the problem solution. 

The correction procedure was repeated until all conditions pertinent to the problem 
were satisfied. Frequent stops were required for adjusting the weighting matrix to avoid 
instabilities in the iteration and to adjust the size of the perturbations so that excursions 
were kept from becoming too large and nonlinear or too small and, hence, inaccurate. 

In general, large weightings were required to hold the initial altitude and velocities close 
to the required values and, hence, avoid divergence of the iteration routine, and small 
weightings were required on the other parameters. 

RESULTS 


is trivial in that the 


quantities must be zero. The inclusion of Xg^tg") 


Numerical Calculations 

Numerical calculations were made to determine the minimum initial mass, two- 
stage rocket vehicle capable of delivering a payload of 4207 kg to a circular orbit about 
the earth at an altitude of 370 400 m (200 int. n. mi.). (See table n for these conditions.) 
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These calculations were made for various 
values of the mass fraction constants as 
shown in table HI to show the changes that 
occur in booster proportions for a range of 
structural weights. Case 1 is a nominal 
design used as a starting point for the cal- 
culations. Cases 1, 2, 3, and 4 corre- 
spond to rather heavy construction. 

Cases 5, 6, and 7 are representative of 
present design practice. Cases 8 and 9 
correspond to very light and nonexistent 
structure, respectively. 


TABLE m.- STRUCTURAL MASS FRACTIONS 
FOR NUMERICAL CALCULATIONS 


Case 

k*p 

k M&e 

k s 

1 

0.10 

0.10 

0 

2 

.10 

.10 

.008 

3 

.08 

.10 

.008 

4 

.08 

.06 

.008 

5 

.06 

.06 

.008 

6 

.06 

.03 

.008 

7 

.03 

.03 

.008 

8 

.03 

0 

.008 

9 

0 

0 

0 


The calculations of the present inves- 
tigation were carried out without regard to 

computational efficiency or sophistication in that a simple differential correction proce- 
dure was employed that made use of numerically determined partial derivatives for iter- 
ation to the point where the necessary conditions were satisfied. Once the first case was 
run and the procedure was established, each succeeding case required about 1 hr of com- 
putation time. As the structural mass fractions were reduced, the problem became more 
sensitive to the choice of the weighting matrix and somewhat more computation time was 
required. One case required in excess of 2 hr of computer time. The adoption of a 
more sophisticated convergence routine including, perhaps, analytically determined par- 
tial derivatives should reduce the machine time. 


All calculations were made for a planar or two-dimensional motion. See the equa- 
tions of appendix B. The rocket was permitted to assume any necessary thrust direction 
for optimization of the vehicle proportions and trajectory. The results of these calcula- 
tions are given in tables IV, V, and VI and in figures 4 to 7. Table IV is a summary of 
the significant results of all of the calculations made. Table V gives the free -space 
velocity increments of the optimized rockets for the various cases. Table VI gives the 
angle of attack of the rocket at maximum dynamic pressure for each case calculated. 
Mass, aerodynamic drag, and velocity are plotted against time in figures 4, 5, and 6, 
respectively, and altitude is plotted against the range angle 9 in figure 7. Case 9 of 
table IQ which corresponds to zero structural mass is not included in these plots. 


Discussion of Results 

Effects on vehicle parameters .- A summary of the results of the numerical calcu- 
lations of optimum two-stage launch vehicles for a 4207 -kg payload and a series of 
structural mass fractions are given in table IV. This table gives the burn time and 
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thrust f oneach of the two stages and, in addition, the initial and final booster and pay- 
load mass. The nominal or starting trial trajectory for these calculations is given 
for reference. 


TABLE IV.- SUMMARY OF RESULTS OF CALCULATIONS OF OPTIMUM LAUNCH VEHICLES 
FOR A 4207-KILOGRAM PAYLOAD DELIVERED TO A CIRCULAR 
EARTH ORBIT AT AN ALTITUDE OF 370 400 METERS 


Case 

Structural mass 
fraction 

Initial 

mass, 

m (ti + ), 

kg 

Mass 

before 

staging, 

Mass 

after 

staging, 

Burnout 

mass, 

“(V). 

kg 

First- 

stage 

burn, 

First- 

stage 

thrust, 

X 7 (tl + ), 

N 

Second- 

stage 

burn, 

Second- 

stage 

thrust, 

x 7( t 2 + )> 

N 


k/p 

k MS e 

k s 

m(t 2 -), 

kg 

m (t 2 + ), 

kg 

*2 ’ tl> 

sec 

*3 - *2’ 

sec 

Nominal 

0.10 

0.10 

0 

646 889 

238 989 

125 463 

23 496 

225.0 

7 116 800 

450.0 

889 600 

1 

.10 

.10 

0 

377 150 

103 908 

32 529 

8 479 

248.8 

4 311 215 

512.6 

184 133 

2 

.10 

.10 

.008 

381 411 

105 572 

33 113 

8 610 

248.1 

4 364 475 

508.5 

189 111 

3 

.08 

.10 

.008 

310 383 

84 659 

31 260 

7 982 

258.1 

3 433 464 

495.8 

184 253 

4 

.08 

.06 

.008 

183 693 

46 510 

22 124 

6 289 

248.8 

2 164 383 

480.3 

129 419 

5 1 

.06 

.06 

.008 

148 507 

45 722 

25 597 

6 377 

179.3 

2 249 767 

462.6 

163 032 

6 

.06 

.03 

.008 

107 177 

27 407 

17 089 

5 268 

176.5 

1 776 393 

435.6 

106 787 

7 

.03 

.03 

.008 

91 649 

25 670 

19 424 

5 034 

192.8 

1 343 064 

459.8 

122 818 

8 

.03 

0 

.008 

73 115 

25 714 

24 094 

4 815 

136.8 

1 359 322 

441.7 

171 345 

9 

0 

0 

0 

62 841 

23 977 

23 977 

4 203 

112.9 

1 350 835 

445.4 

174 259 


An examination of table IV shows a progressive decrease in initial mass of the 
launch vehicle as the structural mass fractions associated with the fuel tanks and rocket 
motors are decreased as was expected. A comparison between the nominal case and the 
rest of the optimized cases indicates that a major effect was the reduction of the thrust 
levels between the nominal and optimized cases, particularly in the second stage. For 
case 1, which has the same structural mass fractions as the nominal case, the second- 
stage thrust was reduced from 889 600 N to 184 133 N. This difference indicates that the 
nominal case was substantially nonoptimal for this particular mission in the sense of 
minimum booster mass. 

A major effect that may be noted as the structural mass fractions are decreased 
in the successive cases calculated (cases 2 to 9) is that the first-stage weight and thrust 
are reduced as the structural mass is reduced but the second -stage thrust and mass 
are held more nearly constant. 

A detailed examination of the results of table IV indicates that the effects on the 
burn times and thrusts resulting from changes in the structural mass fractions are not 
always consistent. For instance, a reduction in k M between cases 3 and 4, as well 
as 5 and 6, causes a reduction in the second-stage thrust Xrj (^ 2 + ) f° r the optimum 
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vehicle; yet a reduction in k^j between cases 7 and 8 causes an increase in the second- 
stage thrust X 7 (t 2 + ). This inconsistency also is evident for the first-stage thrust level 

x 7(*l + )' The burn times t 2 - tj and tg - t 2 are more consistent in respect to this 
change in k M in that they decrease in each case. Some inconsistencies in this same 
regard may be noted for changes in the structural mass factor k T . A complete analy- 
sis has not been made of this point, but some of these inconsistencies may be associated 
with the fact that the vehicle burn time and thrust parameters were not in a sense strongly 
unique near the optimum; that is, changes in the parameters produced little change in the 
initial mass at this point, and continued convergence toward satisfying the necessary con- 
ditions for an optimum became difficult. Such behavior is characteristic of many opti- 
mum solutions, of course, in that the variations of the optimized quantity with respect to 
the state variables approach zero as the optimum is neared. Because of differing sensi- 
tivities and nonlinearities, the necessary conditions of each case calculated were not sat- 
isifed to the same degree. Usually, iteration was continued until little improvement in 
the initial mass could be obtained on successive iterations and until the stability of the 
iteration process became a problem. 

The interstage structural mass fraction kg was changed from 0 to 0.008 between 
cases 1 and 2, and an increase in initial mass was obtained, as would be expected; how- 
ever, the effect of this change was not large. 

Table V gives the velocity increments of the optimized rockets computed on the 
basis of free space (no gravity and no atmosphere) from the appropriate logarithmic 
expressions in mass ratio. The optimized rockets all had a total velocity of between 

TABLE V.- FREE-SPACE VELOCITY INCREMENTS OF 
ROCKETS FOR THE VARIOUS CASES 


Case 

Velocity of 
stage 1, 
m/sec 

Velocity of 
stage 2, 
m/sec 

Total velocity, 
m/sec 

Nominal 

3909 

6576 

10 485 

1 

5060 

5278 

10 339 

2 

5042 

5287 

10 330 

3 

5100 

5358 

10 458 

4 

5392 

4937 

10 329 

5 

4624 

5455 

10 080 

6 

5353 

4619 

9 973 

7 

4995 

5300 

10 296 

8 

4102 

6320 

10 422 

9 

3785 

6831 

10 616 
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9973 and 10 616 m/sec. The optimized rockets, with the exception of cases 8 and 9, all 
hadiapproximately equal velocity increments in each of the two stages. The nominal case 
and cases 8 and 9 had considerably more of the velocity increment in the second stage. 

Effects on trajectories .- Significant trajectory parameters are plotted in figures 4 
to 7 for a number of the cases calculated for this investigation. (See table in.) 

An examination of figure 4 shows a much diminished first-stage mass and fuel- 
consumption rate as the structural mass fractions are decreased. The staging time is 
appreciably lower for the last four cases than for the first three cases shown. A gradual 
consistent trend in the masses and staging times is not evident from these results. Pos- 
sible reasons for this situation have been touched upon in the discussion of the overall 
results of these cases. 


4x10 s 



Figure 4.- Mass, m, as a function of time for cases 2 , 3, 4, 5, 6, 7, and 8. 

An examination of the aerodynamic drag (fig. 5) for these cases shows a similar 
lack of an even or consistent trend as was noted for the mass and staging points in fig- 
ure 4. The maximum drag occurs at various times for the different cases shown, with 
the largest drags occurring for cases 5, 6, and 8. For these cases, this maximum drag 
occurs at an early point in the trajectory where presumably the altitude is somewhat 
low and thus the air density relatively large. The drag impulse or the integral of the 
drag over the time would appear to be more nearly the same for the various cases than 
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the individual details of the drag time history, but even this quantity seems to be consid- 
erably different for the cases shown. 



Figure 5.- Aerodynamic drag as a function of time for cases 2, 3, 4, 5, 6, 7, and 8. 

The velocity time histories for the various cases show a considerable spread and, 
again, give no obvious indication of a consistent trend with structural mass reduction 
(fig. 6). Case 8 shows a larger acceleration rate near the end of its trajectory than do 
the other cases. In this case, the motor structural mass fraction was zero, giving no 
penalty in mass for large thrusts of the rocket motor. This factor may have been signif- 
icant in this particular solution. 

The altitude — range-angle variations for the cases considered show a similar char- 
acter overall (fig. 7). However, case 8 shows a considerably greater loft or steepness 
than the other cases; also, the range angle at orbit injection was considerably less for 
case 8 than for the other cases. This fact seems to be associated, as was noted pre- 
viously, with the greater thrust level and, hence, rate of acceleration for this case than 
for the remainder of the cases. 

An examination of the trajectory printout for case 6 showed that the most severe 
environment insofar as structural bending is concerned occurred at about 70 sec at 
an altitude of 14 420 m, where the aerodynamic angle of attack was 2.61° and the aero- 
dynamic drag was 386 500 N. In this case the product of the dynamic pressure and the 
angle of attack was a maximum; hence, presumably, the aerodynamic bending moment 
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would be a maximum. The aerodynamic angle of 
attack is given by the difference between the flight- 
path angle and the rocket thrust direction (fig. 2). 
The angle of attack at maximum dynamic pressure 
is given in table VI for all cases. For the heavier 
structural masses, cases 1, 2, and 3, the angle of 
attack is as high as 16.9° which is excessive. No 
effort was made to apply a constraint on angle of 
attack for these calculations. 

Many other points of greater or lesser sig- 
nificance could be extracted from the cases calcu- 
lated, but the primary purpose of this investigation 
was the optimization of the significant physical 
parameters of a rocket; therefore, a detailed dis- 
cussion of the trajectories has not been undertaken. 


TABLE VI.- ANGLE OF ATTACK AT 
MAXIMUM DYNAMIC PRESSURE 
FOR THE VARIOUS CASES 


Case 

q max 5 N / m2 

a, deg 

1 

2 176 241 

16.8 

2 

2 203 054 

16.9 

3 

1 830 021 

15.9 

4 

1 571 325 

8.6 

5 

2 407 406 

1.9 

6 

2 446 668 

-.8 

7 

1 787 647 

.2 

8 

2 521 025 

-1.3 

9 

2 952 999 

-2.5 


Extensions to other problems .- A problem such as the optimization of an airplane 
configuration and its flight path may prove to be considerably more difficult to formulate 
than the rocket problem solved in this investigation because the relationships between the 
significant parameters of an airplane and the stage variables and other constants of the 
problem may not be easily determined. No effort has been made to examine an applica- 
tion of this type; however, a few of the more easily established effects dealing with the 
most appropriate or optimum size and sweep of the wings for a given airplane configura- 
tion could likely be determined with a reasonable problem formulation. Such an analysis 
may prove fruitful in that variational methods have not generally been applied to such 
problems. 


CONCLUDING REMARKS 

An analysis using classical variational methods was made whereby both the flight 
path and certain significant physical characteristics of a two-stage rocket were deter- 
mined in an optimal sense. Results were obtained for the most efficient two-stage rocket 
from the standpoint of minimum initial mass that would place a given payload in a circular 
orbit about the earth at an altitude of 370 400 meters (200 int. n. mi.) with a continuous 
burn of the motors. Nine cases were calculated that were representative of a range of 
structural design practices - from very heavy tanks, motors, and interstage structures 
to zero mass structural elements - to show the influence of structural mass. The most 
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^significant result from these numerical calculations was that the second-stage thrust was 
^considerably lower than was expected. Analyses of the type carried out for this investi- 
gation should be useful in the sizing and, generally, in the design of rockets. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., July 31, 1968, 

125-17-05-04-23. 
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APPENDIX A 


EQUATIONS OF MOTION 


The equations of motion of a rocket -powered vehicle operating in the vicinity of the 
earth and in its atmosphere can be derived directly from Newton’s equation of motion 


F = m 


d 2 r 

dt 2 


where the applied forces F are taken to be the gravitational force of the earth, the 
rocket thrust, and the aerodynamic drag, and are given by 


2 = 


^gravitational " m Se ^3 
^rocket = ' Iu = u 


F drag ~ ~ 2 pO*)ACd | f | r 


These quantities are given in relation to the axes shown 
in figure 8, r being the position vector. The drag is 
assumed to be independent of the angle of attack of the 
rocket. No lift force is included in this analysis in that 
the calculations are intended to be illustrative in nature 
and thus have been kept as simple in physical represen- 
tation as possible. 

The equation of motion is then 

m(t)r = u - m(t)g e - |p(h)AC D |r |r 

where m is indicated as a function of time and p is 
a function of altitude. 



The problem is taken to be a planar one that can be adequately represented in polar 
coordinates, and the unit vectors are as shown by 1 and j in figure 8. Resolving this 
equation in terms of these components gives 
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APPENDIX A 


m(t)(r - r0 2 ) = u sin j8 - m(t)g e | p(h)AC D [i 2 + (re) 2 ] ^ r 

m(t)(re + 2re) = u cos /3 - i p(h)AC D [l 2 + (re) 2 ] ^ re 




(Al) 


where /3 refers to the inclination of the thrust to the local horizontal (fig. 2). In this 
problem an exponential representation of the earth atmosphere is used as given by 

-Wr-r ) 

p(h) = p(r - r e ) = p e e ' ' 

rather than a tabular form in order to simplify the determination of p. In this repre- 
sentation v is a constant and p e is the mass density of the atmosphere at the earth 
surface (ref. 7). The flight -path angle y is given as 

y = tan~l — 

rd 
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APPENDIX B 


THE NECESSARY MULTIPLIER RULE 

In reference 2, a multiplier rule is formulated for a Bolza type variational problem 
that was later extended in reference 3. This multiplier rule as extended may be applied 
to the present problem to derive necessary conditions for a minimizing arc. This rule 
corresponds to the multiplier rule studied by Bliss (ref. 1) and others except that the 
function to be minimized depends not only on the coordinates of the end points of an 
admissible curve, but also on the coordinates and discontinuities at variable intermediate 
points on this curve. The coordinates of the end points and the discontinuities in coordi- 
nates or the coordinates at intermediate points are further restricted to satisfy certain 
side conditions. A short description of this rule is given here for completeness along 
with a brief comment on the use of the necessary condition on the Weierstrass E -function. 

This multiplier rule concerns itself with finding in a class of sets Jt a , x(t)J which 
satisfies the differential equations 

0. [t,x,x] = 0 (j = 1, . . ., m < n) 


and the end and corner conditions 

£a> x i( t a + )> x i( t a‘)] = 0 (p = 1, . . ., r ^ (n + 1) (q + 1)) 

(where the superscripts + and - indicate conditions immediately after and prior to the 
designated point in time), a set which minimizes a functional J of the form 

J = g[ta,x i (t a + ),x i (t a -)] 

where t a denotes the set of variables t 2 , . . ., tq) which is always understood to 
satisfy the inequalities tj < t£ < . . . tq and Xj(t), or simply x(t), denotes the set of 
functions 

[^(t), x 2 (t), . . ., x n+1 (t)] 

In these expressions, Xj(t a ") and Xj(t a + ) denote the values of x^ 

each side of t a where a discontinuity or corner may exist, and Xi(t) 

(i = 1, . . ., n + 1). For the present problem, t a will be taken as the 


(*1 5 ‘ = ‘ q ) 

immediately on 
_ dXi(t) 
dt 
time. 
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APPENDIX B 


Anradmissible solution E of the equations = 0, defined on an interval |t^,tqj, 
is saidrto: satisfy the multiplier rule if there exists a function 

G [^a^il^-a )> x i (^a + )’^0^/iJ = ^0^ fa ,x i (*a )j x i(V + )J + {*a ,x i(*a )’ x i(^a + )J 

with constant Xg, not all equal to zero, and a function 

F(t,x,x,X) = Xj(t)0j 

with multipliers Xj(t) continuous on t^,tq except possibly at corners of E where 
they have well-defined right and left limits, such that the following equations are 
satisfied: 



where the subscript c designates a corner or point of discontinuity; a number of these 
may exist. For an arc E satisfying the multiplier rule, the multipliers Xq and X^(t) 
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APPENDIX B 


do not vanish simultaneously at any set of associated points on E. Every minimizing 
arc must satisfy the multiplier rule. 

The necessary condition on the E function as given by Bliss (ref. 1), which is 
equivalent to requiring that the function H = be a maximum (ref. 8), may be used to 

choose the minimizing arc from among those satisfying the multiplier rule. 
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